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A number of coupled particle- element and hybrid particle- element methods have been de- 
veloped for the simulation of hypervelocity impact problems, to avoid certain disadvantages 
associated with the use of pure continuum based or pure particle based methods. To date 
these methods have employed spherical particles. In recent work a hybrid formulation has 
been extended to the ellipsoidal particle case. A model formulation approach based on La- 
grange’s equations, with particles entropies serving as generalized coordinates, avoids the 
angular momentum conservation problems which have been reported with ellipsoidal smooth 
particle hydrodynamics models. 
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INTRODUCTION 

A review of the literature on hypervelocity impact simulation shows that most research in 
this field has focused on the development and application of continuum hydrocodes, using ei- 
ther Eulerian, Lagrangian, or Arbitrary Lagrangian-Eulerian (ALE) formulations [1,2,3]. In 
the last decade attention has shifted towards particle based methods [4], in particular smooth 
particle hydrodynamics (SPH). Although both continuum and particle based methods have 
demonstrated excellent results in a range of practical applications, these methods are not 
without problems. As a result some recent research has formulated coupled particle-element 
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[5] or hybrid particle-element [6] methods, aimed at avoiding certain disadvantages of the use 
of pure continuum based methods or pure particle based methods in hypervelocity impact ap- 
plications [7]. The present paper develops a generalized hybrid method, extending the work 
of Fahrenthold and Horban [6] , presenting for the first time a combined particle-element for- 
mulation based on a nonspherical particle geometry. Nonspherical particle formulations are 
also of interest in certain molecular dynamics [8] and astrophysics [9,10] applications, where 
material or flowheld geometry has suggested the use of repulsion potentials or interpolation 
kernels with ellipsoidal forms. 

Most previous hypervelocity impact simulation work has employed Eulerian finite volume 
and Lagrangian finite element methods [1,11] and more recent extensions of these methods 
to Arbitrary Lagrangian- Eulerian frames [12]. Although continuum hydrocodes are accurate 
and efficient simulation tools for many important problems, they are not well suited for use in 
all hypervelocity impact applications. Eulerian codes can accommodate arbitrarily general 
contact-impact, but their approximate models of material strength effects are best suited to 
a very high velocity impact regime. Lagrangian codes incorporate very accurate models of 
material strength effects, but their slideline based contact-impact algorithms are best suited 
to a relatively low velocity impact regime. ALE methods allow for intelligent tailoring of the 
mesh, ideally avoiding (but not eliminating) the preceding limitations. 

None of the aforementioned continuum formulations is ideally suited to applications where 
the generation, transport, and contact-impact dynamics of material fragments is of central 
concern. When material fragments become much smaller than the finite volume cell size, 
Eulerian flow held interpolations may not accurately represent the physical state of the highly 
comminuted medium. Hence the adaptive introduction of a high resolution Eulerian mesh 
is required. Similarly the use of Lagrangian finite elements to represent a highly fragmented 
medium calls for the introduction of a very extensive and adaptive slideline mesh. It is 
not surprising that these continuum formulations can be difficult to apply in cases where 
fragmentation dynamics are of central interest, given the highly discontinuous nature of the 
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field variables in such systems. 

Recognizing the disadvantages of continuum formulations in certain hypervelocity im- 
pact applications, a number of different particle methods [13,14,15] have been introduced. 
In general these methods replace or augment the conventional continuum kinematics with 
particle based interpolations, to more efficiently represent fragment generation, transport, 
and contact-impact effects in applications where modeling of the latter physics is of central 
concern. The overwhelming majority of particle based hypervelocity impact models have 
employed SPH techniques and spherical kernels. A notable exception is the work of Shapiro 
and co-workers [9,10], who extended the basic SPH method to ellipsoidal kernels. However 
those authors indicate that their method fails to conserve angular momentum, and they of- 
fer no solutions to the tensile instability, numerical fracture, and other problems which have 
hindered the effective use of SPH methods in some applications. The formulation presented 
here does not involve element-to-particle conversion, as discussed for example by Johnson et 
al. [14], who converted distorted finite elements into SPH particles. The latter conversion, 
if it occurs before failure of the element, may introduce the aforementioned SPH tensile 
instability and and numerical fracture problems. 

The present paper formulates a hybrid ellipsoidal particle-finite element model for hyper- 
velocity impact simulation. Unlike the coupled particle-element formulations of Hayhurst et 
al. [16] and other authors, which use particles to model some structures and finite elements 
to model other structures, the present paper employs particles and elements everywhere. The 
simultaneous use of particles and elements is not redundant, since they are used to represent 
distinct physics. The particles model all inertia and contact-impact effects, and are associ- 
ated with interpolation functions which represent compressed states. The elements model all 
strength effects, namely tension and elastic-plastic shear. This approach has several advan- 
tages. The tensile instability and numerical fracture problems common to SPH formulations 
are avoided, since interpolation kernels are not used to represent interparticle tension forces, 
and since material strength effects give rise to interparticle forces only among reference con- 
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figuration neighbors. Although true Lagrangian kinematics are used to calculate tensile and 
shear forces, the use of particles in all inertia and contact-impact calculations means that 
material failure is simply accommodated (without element-to-particle remapping or mass 
and energy discard) via the loss of element cohesion. Particles not associated with any in- 
tact elements translate and rotate as individual fragments, in response to contact-impact 
loads. No slidelines are used, so that contact-impact of all intact and fragmented material 
is modeled everywhere. Although the authors do not contend that this numerical approach 
is best for all impact simulation problems, it offers important advantages in hypervelocity 
impact applications. In the latter field an interest in multi-structure perforation, erosion, 
fragmentation, and very general contact-impact modeling is not unusual. 

The model formulation work described in the sections which follow is different from that 
normally employed to develop hydrocode models of either the continuum or particle type. 
Two particular features should be mentioned. First an energy method (Lagrange’s equations) 
is used to develop the semi-discrete model, so that no reference is made to any partial 
differential equations. The introduction of entropy variables as generalized coordinates makes 
it possible to account for very general thermal dynamics while using a discrete Lagrangian 
approach. Secondly the formulation introduces Euler parameter states [17] to account for 
the rotational motion of the nonspherical particles. An Euler parameter description of the 
particle kinematics is desired in order to avoid the singularities of Euler angle based models. 
Although the use of Euler parameters leads nominally to a differential-algebraic formulation, 
the development which follows introduces momentum state variables, to reduce the order 
of the system and thereby eliminate the Lagrange multipliers associated with the Euler 
parameter constraint. The model formulation work described here therefore generalizes in 
a geometric sense the particle-element work of Fahrenthold and Horban [6] and generalizes 
in a thermomechanical sense the Euler parameter based mechanical models of Chang and 
Chou [18] and others [19,20]. 

The organization of this paper may be summarized as follows. The second and third 
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sections define the particle and element kinematics and the interpolations used for all held 
variables. The fourth and fifth sections develop kinetic energy and thermomechanical po- 
tential energy functions for the system. The sixth section describes the dissipative process 
models, for plastic formation and damage evolution. The seventh and eighth sections present 
artificial viscosity and artificial heat diffusion models, similar to those used in other shock 
physics codes, and develop entropy evolution equations which take the place of the energy 
conservation relations normally employed. The ninth section introduces a virtual work ex- 
pression, to account for externally applied loads. The tenth section combines the canonical 
thermomechanical Lagrange equations with the multiple nonholonomic constraints devel- 
oped in the preceding sections, and derives an unconstrained explicit state space model 
for the particle-element system. The final section presents two example problems, show- 
ing good agreement of simulations performed using the model developed here to published 
data for three dimensional impact experiments. Additional validation work, for simulations 
performed using a spherical kernel, are described by Fahrenthold and Shivarama [21], who 
also provide details on the numerical implementation and test results showing numerical 
convergence and good parallel speedup. 


PARTICLE KINEMATICS 


The inertia distribution in the modeled system is represented by a collection of n ellipsoidal 
particles, with the mass of the Ah particle, whose position and orientation are described 
by a center of mass position vector (c^O) and an Euler parameter vector (e^) with component 
forms 


(0 - r Jd AO AO i t 


c = c 


(0 = [ e (0 ao AO AO } T 


e K ' = 


o c i 


( 1 ) 


where T denotes the transpose. The components of the center of mass position vector 
are described in a fixed Cartesian coordinate system. The Euler parameters [17] define an 
orthogonal rotation matrix (R(0) for each particle 


R« = hW g wt 
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which relates the components of any vector (a), described in a body fixed co-rotating coor- 
dinate system, to its corresponding components in the fixed Cartesian system, using 

a = R® a (5) 


The four Euler parameters may be used to compute three Euler angles for each particle, 
and are preferred for use in large rotation dynamics problems, since the Euler angles are 
singular functions. Although the Euler parameters are nonsingular, they must satisfy the 
constraint 

e® T e® = 1 ( 6 ) 

which requires that 

G®e® = 0, G®e® = -G®e®, G® G® T = ! (7) 


where I denotes an identity matrix. 

The Euler parameters and their time derivatives are related to the particle angular velocity 
vectors ( u>® ), with components expressed in the particle’s co-rotating frame, by the well 
known [17] relations 

= 2 G® e®, e® = ^ G® T u>® (8) 

Similarly the anti-symmetric matrix ft, with axial vector u>, which satisfies 

G®v = w (0 x v ( 9 ) 


for any vector v, is related to the Euler parameters and their time derivatives by 

ft ® = 2 G® G® r = -2 G®G® T (10) 
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The next section presents a density interpolation associated with the particles just de- 
scribed, and defines Lagrangian finite elements whose nodal coordinates are components of 
the particle center of mass position vectors. 


DENSITY INTERPOLATION AND FINITE ELEMENTS 


In the material reference configuration, the ellipsoidal particles are arranged in a packing 
scheme defined by a mapping from a body centered cubic unit cell, the mapping composed 
of stretchings in three orthogonal directions, those directions aligned with the principal axes 
of the ellipsoids. The density interpolation for compressed states is 


p ( i) = p f + ^ iy«> (ii) 

j = 1 


where is the continuum density at the centroid of the fth particle, pd ! is the reference 
density for the ith particle, and is the interpolation kernel 
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The symbols S tJ and u denote respectively the Dirac delta function and the unit step function 


u(x ) 


0, x < 0 

1, x > 0 


(14) 


The constants a and (3 are effective separation distances, in a unit cell, measured respectively 
between body centered particles and between a body centered particle and a particle located 
at a cell vertex 


a = 




(15) 


Note that the lead coefficient in the expression for is due to the presence of eight 

nearest neighbors in the reference configuration. 


The function is an ellipsoidal coordinate, defined in a frame which co-rotates with 


the jth particle 
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where h\ 3 \ h^\ h\p are the half-lengths of the principal axes of the j'th particle. 

Several properties of the density interpolation and kernel functions just described should 
be noted. The interpolation provides an exact Lagrangian description of the variation of 
density under uniform compression, and incorporates a lower bound equal to the particle ref- 
erence density, so that tensile instabilities are avoided. The kernel is singular, so that no two 
particle centers overlap, avoiding particle streaming effects. In the reference configuration 
neighbor particles do not contribute to the density calculation, so that no special treatment 
of surface particles with incomplete neighbor sets is required. The density dependence of 
the compact kernel support means that particles interact only with near neighbors, and that 
mechanical interaction with remote neighbors through intervening matter is avoided. Finally 
note that the dimensionless kernel used here is closer in functional form and physical inter- 
pretation to the repulsion potentials used in molecular dynamics, than to the dimensional 
kernels with cubic spline form used in most SPH formulations. 

The center of mass coordinates of vertex centered particles are nodal coordinates for hex- 
ahedral finite elements, which are used to account for tension and elastic-plastic shear in 
cohesive solid materials. These elements are used with one point integration, as described 
by Hallquist [22], although in the present case the element packing scheme and density in- 
terpolation preclude the development of hourglass modes. Since large deformation problems 
are of interest in hypervelocity impact applications, the element nodal coordinates are used 
here to compute a Jacobian (J^) and a Lagrangian deviatoric strain tensor (E^) for each 
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element, defined by [23] 


Jti) = det F (j) , E U) = ^ ( C (i) - 1 ) (18) 

where F (y ^ is the deformation gradient for the jth element and 

C°) = F°' )T F 0) , F {j) = ( det F F W) (19) 

The plasticity model discussed in a later section assumes an additive decomposition of the 
deviatoric strain, with the elastic deviatoric strain tensor (E ft W) f or the jth element dehned 

by 

E e(j) = E 0) - E P(J ') (20) 

so that the evolution equations for the plastic strain tensor (E p bl) must satisfy the isochoric 
plastic deformation constraint 

tr (c p( b- T C p(j) ^ = 0, E p(i ) = I (cpCj): _ i) (21) 

It should be emphasized that the numerical modeling methodology developed in this paper 
applies for a wide range of element types and plasticity models. The particular element 
kinematics and constitutive equations used in this paper are representative formulations 
which account for large deformation kinematics. 

KINETIC ENERGY 

The kinetic co-energy for the system is 

n 

T* = ^2 T * (i) ( 22 ) 

1 = 1 

where T*W j s the kinetic co-energy for the ith particle 

T*(*) = - m {i) c w T c w + - T jW w (i) (23) 

2 2 v ’ 

and the components of the inertia matrix Jb) are calculated in a body fixed frame. Since 

the components of the angular velocity vector are quasi-velocities, equations (8) are used to 
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rewrite the second term 




- c ^ T c« 
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+ 2e (i)T G (i)T J (i) G w e w 


(24) 


which identifies the center of mass coordinates and the Enler parameters as generalized 
coordinates in Lagrange’s equations. Note that the relations (7) allow the rotational kinetic 
co-energy to be expressed as 


rp*(i ) 
rot 


2e (i)T G (i) T J (i) G w e w 


(25) 


so that the generalized force associated with the Euler parameter dependence of the kinetic 
co-energy is 

r)T*(h 

k h) = = 4 G« t J«g« e (i) (26) 

<9eW v ’ 

The next section develops a potential energy expression for the particle-element system. 


POTENTIAL ENERGY 

The potential energy of the particle-element system is the sum of the particle internal 
energies and the element stored energies due to tension and elastic shear, and has the general 
form 


^ m(i) u{i) + J 2 ( £/ten(j) jjdev(j) _|_ jjpen(j) j 

3 = 1 


(27) 


i= 1 


where j s an internal energy per unit mass, n e is the number of elements, U ten ^ is a 
stored energy in tension, JJ dev ^' ) is a stored energy in shear, and U pen ^ is a penalty energy 
used to position the body centered particle in each element. The particle internal energy 
density is determined by the mass density and entropy density (s^) 


U® =U {i) (p W ,S W ) , 


s « = 


m 


(0 


(28) 


where S ^ is the total entropy of the fth particle. The functional form of the internal energy 
per unit mass depends on the equation of state, and it is not unusual in such calculations to 
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rely upon tabular data [24]. In any case, the extensive state variables which determine the 
particle internal energy are 

u® =u® (c®,e®,S®) (29) 

The form of the element stored energy functions depends upon the constitutive assump- 
tions. Although the general method developed here admits a wide range of constitutive 
models, simple forms are assumed here, and are described as follows. The stored energy in 
tension is 

U ten U) = I (l - D u) ) k& V} j) (J u) - l) 2 u{ J u) - 1) (30) 

where is a bulk modulus, is the normal damage, and vd 3> is the element reference 
volume. The stored energy in shear is 

U dev{j) = (1 - d ij) ) fi ij) Vj j) tr (E e ^ T E e(j ' } ) (31) 

where ii ( j> is a shear modulus and d {j) is the deviatoric damage. The penalty function is 

U pen{j) = ^ (1 - D( j) ) K® (c (j) - c (i) ) 2 (32) 

where is the center of mass position vector for the body centered particle of the jth 
element, is the average of the center of mass position vectors for the vertex centered 
particles (the element nodes), and the penalty stiffness is 

K {j) = E U) V 0 (j) ^ (33) 

where is Young’s modulus for the element. The damage variables are used to model the 
loss of cohesion due to element failure, and are discussed in the next section. 

The preceding constitutive assumptions combine with the element kinematic relations 

jU) = jU) ( C M) ? E e(j) = E e(j) epOO) (34) 

to yield the following state variable dependence for the system potential energy 

V = V (c w , e w , S (i) , E p(j \ D®, d U) ) (35) 
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The listed arguments, which include the particle entropies, are by definition Lagrangian 
generalized coordinates. 

The system potential energy function defines the conservative generalized forces 


dV 


= g« 


dV 
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dc® ’ <9eb) ’ dS® 

where 6® is a particle temperature, and gW and are forces and torques which act on 
the ellipsoidal particles. The energy conjugates for the dissipative state variables define a 
deviatoric stress 


as well as the energy release rates 

Y d U) = 
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The next section discusses evolution equations for the dissipative state variables. 


(37) 


(38) 


PLASTICITY AND DAMAGE MODELS 

The evolution equations for the plastic strain tensor and continuum damage variables will 
serve as nonholonomic constraints on the system level Lagrange equations. As discussed in 
the last section, the general formulation developed here admits a wide range of constitutive 
models. The relatively simple dissipative constitutive models assumed here are described 
in this section. The plasticity model is adapted from Fahrenthold and Horban [25], and 
represents the simplest possible accommodation of the isochoric deformation constraint (21). 
The damage evolution relations are adapted from the Eulerian hydrocode work of Silling [26]. 

The flow rule for the plastic strain is 


j. 0) 
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where A^) is a scalar multiplier and 
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(39) 

(40) 

(41) 


13 


e pU) = 


i tr (^pU)T-^p(j)^ 


1/2 


( 42 ) 


with rj (j> a strain hardening modulus, a strain hardening exponent, and e p ^ the accu- 
mulated plastic strain. The effective stress tensor which appears in the flow rule is 


gpd) _ ^ _|_ ^d) € pU)\- nW gd) ; 


s “ = ^^) = (1 - d “ )2 ' , “ E, “ (43) 
where the indicated deviatoric stress is power density conjugate to the plastic strain rate, in 
the entropy equality for the solid medium. 

The yield condition is 


yd) _ T d) _ yd) 


T d) _ 


-tr (S p(i)T S p(i) ) 


1 1/2 


( 44 ) 


where r ^' 1 is the second invariant of the effective stress. The yield stress Y is 

y(i) = yd) (1 _ 7 (J)0«(J)) 


( 45 ) 


where is the reference yield stress, 7^ is a thermal softening modulus, and 0 H ^ is the 
maximum historical homologous temperature. The plastic strain increment at each time 
step is determined using a one step iteration procedure [22] with 

-O') _ yd)) fi( T d) _ yd)) 


AA (i) = 


T v 


(1 - dd)) 2 pd) 

The final evolution equations for the plastic strain have the functional form 


( 46 ) 


E rU) _ j;rU) 


( 47 ) 


and are nonholonomic constraints on the particle-element model. 

Hypervelocity impact problems normally involve perforation and fragmentation effects, so 
that the ability to model such physics is essential. I11 the present case normal and deviatoric 
damage variables are introduced, to model the transition from a cohesive solid, characterized 
by intact finite elements, to a comminuted medium, described by the free flow of particles. 
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Free particles interact with cohesive material and with other fragments under general contact- 
impact loads. 

The damage evolutions applied here are 

D® = u( 1 - D (j) ), d ® u( 1 - d (j) ) (48) 

fiAt K h n At v v ' 

where At is the time step, n is the number of time steps used to model the transition from 

an intact to a failed element, and 

A w = max{ u(e p{3) - e p f U) ), u(6^ ax - 9$), u(J [ c j) - J^ n ), u(a & L - } (49) 

The function just defined initiates damage evolution in an element when when the accumu- 
lated plastic strain e p( ^\ maximum temperature 9 max , minimum element Jacobian J ^] n , or 
maximum eigenvalue of the deviatoric stress (Jrnax reach corresponding critical values for the 
failure strain (e^) , melt temprature (6m), maximum compression (jj^), or spall stress 
(aF' ) ). More sophisticated damage evolution relations may be introduced without change 
to the general modeling framework. The damage evolution equations are nonholonomic 
constraints which apply to the system level Lagrange equations. 


ARTIFICIAL VISCOSITY AND HEAT DIFFUSION 


The shock physics problems of interest here call for the introduction of artificial viscosity 
and artificial heat diffusion effects. The standard formulation used here takes the viscosity 
force pb) on the ith particle to be 


pw = ^2 u(i,j) max (°y ij) ) r M , 

3 = i 


r (*d) 
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where the velocity difference ybJ) j s positive for converging particles 
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where ci l> is the sound speed for the ith particle and the parameters c Q and c\ are dimen- 
sionless linear and quadratic numerical viscosity coefficients. 

The numerical heat diffusion model assumed here takes the thermal power flow for the 
ith particle to be 

n 




1 

2 


QCon(i) 


Y R (i,j) ( 0 (<) ~ 0 {j) ) 

3 = 1 



(53) 

(54) 


where is the specific heat for the it\i particle and the parameter k Q is a dimensionless 
numerical heat diffusion coefficient. 


ENTROPY EVOLUTION EQUATIONS 

The introduction of entropy state variables, which facilitates the use of an energy based 
modeling approach, means that the conservation of energy relations normally incorporated 
in shock physics models are replaced here by evolution equations for the particle entropies. 
These evolution equations are 

i ) nirr(i) 


gcon(i) 


(55) 


where the irreversible entropy production (S'* rr M) is due to plastic deformation, damage 
evolution, and viscous disspation 




f»( i ) T c(*) 


+ E^‘ 


<j) 


+ r + V„ w tr (s u)T E pW 


(56) 


3 = 1 J 

and is the fraction of the dissipation in element j associated with the ith particle. The 
conduction entropy flow ( 1 S' con (*)) is due to artificial heat diffusion 

n gti) 

gcon(i) = 9 (i)-l QCon(i) = ^ ( 1 - — ) (57) 

3 = 1 

Since the dissipated energy is not lost but rather transduced to the thermal domian, the en- 
tropy evolution equations are nonholonomic constraints in the thermomechanical Lagrangian 
formulation. 
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VIRTUAL WORK 


Although external loads are normally not considered in hypervelocity impact applications, 
for completeness this section develops a virtual work expression. The virtual work for the 
system is a summation over the particles 

n 

sw = 6w(i) (58) 

i = 1 

In this case the particle angular velocities are the time derivatives of rotational quasi- 
coordinates (qW) 

q W = (59) 

so that the virtual work for the ith particle is 

5W {i) = f (i)T hc w + T (t)T hq w (60) 

where and TO are externally applied forces and torques. In view of equation (8), virtual 
changes in the quasi-coordinates and the Euler parameters are related by 

hq (i) = 2 G w he w (61) 

so that the particle virtual work expression is 

5W (i) = f 1 ' i)T 5c {i) + 2 [ G wt T w ] T he w (62) 

The coefficients of the virtual changes in the generalized coordinates are by definition gen- 
eralized nonconservative forces in the system level Lagrange equations. 

LAGRANGE’S EQUATIONS 

The stored energy functions, constraint equations, and virtual work expression developed 
in the preceding sections may be combined with the canonical Lagrange equations, to obtain 
an ODE model for the particle-element system. 


17 


The canonical Lagrange equations are 


d / dT* \ dT* 
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where Q c ^, Q e ^, Q s ®, Q d ®, Q D ^\ and Q p ® are generalized forces determined by the 
constraints and the virtual work. Note that the rate form of the Euler parameter constraint 
for the it\i particle is 

e (i)T e (i) = 0 (67) 

Introducing a Lagrange multiplier j e ® for each Euler parameter constraint as well as La- 
grange multipliers 'y d ® j ^d{i) ^ anc [ yvd) f or the entropy, normal damage, deviatoric 
damage, and plastic strain evolutions equations, the generalized nonconservative forces are 
found to be 
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If these results and the degenerate Lagrange equations (65) and (66) are combined with the 
previously derived expression for the kinetic co-energy, the momentum balance relations take 
the form 


d 

dt 


( m ® c^) + 


dV 

<9cb) 


f (*) _ p(») 


(74) 


d 

dt 


(4G (i)r J w G w e (i) ) -k w + 


dV 

<9eh) 


2 QbA'jM ^e(i) e (i) 


Introducing the momentum variables 


p ('0 = m (0 c w 


h w = J w 


(75) 


(76) 


and premultiplying the angular momentum balance expression (75) by | G^, which [using 
equation (7)] removes the last term, yields the final state space formulation 


pW 

_ _g(») _ p(») _|_ f(<) 

(77) 

hW 

= -f2 (i) h (i) - -G (i) M (i) + T (i) 
2 

(78) 

c« 

= m (<)_1 p (<) 

(79) 

e« 

— -G( i ) T j(*) -1 h( i ) 
2 

(80) 


Augmented by the evolution equations for S^'\ D (j \ d! J \ and E p b), these are explicit first 
order equations for the particle-element system. With appropriate initial conditions they 
may be integrated to obtain solutions for a wide range of hypervelocity impact problems. 


EXAMPLE SIMULATIONS 

This section compares the results of simulations performed using a numerical implementa- 
tion [27] of the model developed here, to published experimental results for two hypervelocity 
impact problems. The simulations employed a Mie-Grunsisen equation of state and 

c 0 = 0.001, ci = 0.0, k Q = 0.1, e p f = 1.0 (81) 
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with material properties taken from Steinberg [28]. 

The first example involves the impact of a uranium alloy rod on a steel plate, 0.64 cm in 
thickness, in an experiment described by Hertel [29]. The cylindrical projectile is 0.767 cm 
in diameter, with a length to diameter ratio of ten, and impacts the plate at an obliquity of 
73.5 degrees and a velocity of 1.21 kilometers per second. The velocity of the target plate is 
0.217 kilometers per second. The simulation employed 712,929 particles, with the projectile 
composed of spherical particles and the target plate composed of ellipsoidal particles with 
aspect ratios of 1.5:1. 5:1.0. Figure 1 shows an element plot of the initial configuration, while 
Figure 2 shows a particle plot of the simulation results at 100 microseconds after impact. 
The simulation results for the rod erosion (27.5 percent) and residual rod velocity (1.10 
km/s) show good agreement with the corresponding experimental values (27.6 percent and 
1.07 km/s). 

The second example involves the impact of a tungsten rod on a steel plate, 0.95 cm in 
thickness, in an experiment described by Yatteau et al. [30]. The cylindrical projectile is 
0.475 cm in diameter, with a length to diameter ratio of twenty, and impacts the plate at 
an obliquity of 75 degrees and a velocity of 1.83 kilometers per second. The target plate 
is stationary. The simulation employed 671,176 particles, with the projectile composed of 
spherical particles and the target plate composed of ellipsoidal particles with aspect ratios of 
1.5:1. 5:1.0. Figure 3 shows an element plot of the initial configuration, while Figure 4 shows 
a particle plot of the simulation results at 150 microseconds after impact. The simulation 
results for the rod erosion (37.9 percent) and residual rod velocity (1.60 km/s) show good 
agreement with the corresponding experimental values (40 percent and 1.78 km/s). 

In the preceding problems the vast majority of the particles are associated with the target 
plates, so that the use of ellipsoidal particles reduces the total particle count by a factor of 
approximately 2.25, equal to the product of the three indicated aspect ratios. It should be 
noted that the computational advantages of a reduced particle count can in part be offset 
by an increase in the average number of nominal neighbors per particle. This increase is 
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due to the fact that the neighbor search for each particle must consider a spatial volume 
proportional to the cube of the major ellipsoidal axis, despite the fact that a closer approach 
without contact can occur when particle minor axes are appropriately aligned. Experience to 
date indicates that the use of ellipsoidal particles leads to a reduction in memory requirements 
which is proportional to the reduction in particle count, but to a cpu time requirement similar 
to that measured in corresponding simulations with spherical particles. 

CONCLUSION 

The present paper has developed an ellipsoidal particle-finite element method for hyper- 
velocity impact problems, and demonstrated its use in the simulation of three dimensional 
impact experiments. The hybrid particle-element kinematic scheme allows for the use of 
true Lagrangian strength models, while incorporating a completely general description of 
contact-impact effects. Tensile instability, numerical fracture, and angular momentum bal- 
ance problems, often associated with particle formulations, are avoided. Material remaps, 
mass or energy discard, slideline algorithms, and adaptive mesh refinement, often associated 
with continuum models, are also avoided. Derivation of the model is based on a thermome- 
chanical form of Lagrange’s equations, with entropy variables as generalized coordinates. The 
use of an energy based modeling approach facilitates the systematic integration of diverse 
particle based and element based interpolations. Work to date suggests that the method pro- 
vides a numerically robust and accurate approach to the simulation of hypervelocity impact 
problems. 
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Figure 1. Uranium alloy long rod impact on a steel plate, element plot of the initial 
configuration 

Figure 2. Uranium alloy long rod impact on a steel plate, particle plot at 100 microsec- 
onds after impact 

Figure 3. Tungsten long rod impact on a steel plate, element plot of the initial config- 
uration 

Figure 4. Tungsten long rod impact on a steel plate, particle plot at 150 microseconds 
after impact 
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Figure 1: Uranium alloy long rod impact on a steel plate, element plot of the initial config- 
uration 
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Figure 2: Uranium alloy long rod impact on a steel plate, particle plot at 100 microseconds 
after impact 
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Figure 3: Tungsten long rod impact on a steel plate, element plot of the initial configuration 
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Figure 4: Tungsten long rod impact on a steel plate, particle plot at 150 microseconds after 
impact 
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